Diffusion and criticality in undoped graphene with resonant scatterers 
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A general theory is developed to describe graphene with arbitrary number of isolated impurities. 
The theory provides a basis for an efficient numerical analysis of the charge transport and is applied 
to calculate the minimal conductivity, o, of graphene with resonant scatterers. In the case of smooth 
resonant impurities a grows logarithmically with increasing impurity concentration, in agreement 
with renormalization group analysis for the symmetry class DHL For vacancies (or strong on-site 
potential impurities) o saturates at a constant value that depends on the vacancy distribution among 
two sublattices as expected for the symmetry class BDI. 
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Transport properties of graphene [IMS'] remain in the 
focus of intense studies. It has been established, both 
theoretically [1, Q and experimentally that the 

conductivity of short and wide samples of ballistic 
graphene acquires a minimal value of 4 x /irh (the fac- 
tor 4 reflecting the spin and valley degeneracy) when the 
chemical potential is tuned to a vicinity of the Dirac 
point. The minimal conductivity of larger graphene 
flakes is close to 4 x e^//i for the majority of experimen- 
tally available samples [2] . The enhancement of the min- 
imal conductivity is attributed to the effect of disorder: 
graphene near the Dirac point may conduct better when 
impurities are added lll-[l3{. 



Remarkably, the minimal conductivity of disordered 
graphene remains essentially constant when the temper- 
ature T is lowered by several orders of magnitude, down 
14| . This is in contrast with the behavior of 

32//^: 



to 30 mK 

conventional 2D systems with conductivity a 



their a gets strongly suppressed with lowering T due 
to Anderson localization. The absence of localization in 
graphene indicates that the dominant disorder is cither 
of long-range character (and thus does not mix the val- 
leys )or preserves a chiral symmetry of the Hamiltonian 



^ .^The former possibility has been investigated in 

Refs. |l2. 13, 16, 17 1 . In this paper we explore the case 
of resonant impurities that preserve the chiral symmetry 
[Cz in terminology of Ref. 15]). 

Resonant scatterers create "mid-gap" states directly at 
the Dirac point, thus having a strong impact on the min- 
imal conductivity. A natural example of a resonant scat- 
terer is a strong potential applied to a site of a graphene 
honeycomb lattice, which is equivalent to a vacancy. In 
this way, vacancies are effectively created by hydrogen 
ad-atoms or CII3 molecules, which bind to a single car- 
bon atom in graphene and change its hybridization from 
sp^ to sp^ type. The resonant character of hydrogen 



ad-atoms was supported by the DFT analysis jl8|. Res- 
onant scatterers give linear (up to a logarithmic factor) 
dependence of the graphene conductivity on electron den- 
sity 1^ l^j consistent with experimental observations. 
Recent experiments (2^, [2l| confirmed that a moderate 
concentration of hydrogen adsorbates preserves all the 
salient features of transport in graphene and provided ev- 
idence that resonant impurities determine the graphene 
mobility. 

The vacancies are not the only type of resonant impu- 
rities. For instance, a smooth potential impurity, which 
can be represented by a scalar potential in the Dirac 
equation, is resonant provided the energy of the local- 
ized impurity state coincides with the Dirac point. Re- 
cent works 22, 23] studied the effect of resonant scalar 
impurities in the ballistic regime. It was shown, in partic- 
ular, that each such impurity enhances the conductance 
of ballistic graphene by a value of the order of /h. 

In this Letter we develop a general theory of trans- 
port in a system with arbitrary number N of isolated 
impurities. The full counting statistics of the system 
is given by a determinant of a matrix of size N . The 
averaging over impurity positions can be performed nu- 
merically with a high efficiency. We apply this theory 
to disordered graphene at the Dirac point, focusing on 
two types of resonant scatterers: scalar (smooth poten- 
tial) impurities and vacancies. This allows us to explore 
the diffusive (or critical) regime that is established with 
increasing concentration of impurities (or, equivalently, 
sample length). 

We consider a graphene sample of the length L and the 
width W ^ L, which is described by the Dirac Hamil- 
tonian, H = —ihvaV + V{y), where cr = {ax, cry) is the 
vector of Pauli matrices, v is the velocity, and V{r) is an 
impurity potential that can mix valleys and sublattices. 
Hereafter we set Hv — 1. 
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Metallic leads at a; < and x > L are defined by 
adding a large chemical potential, /loo oo. Inside the 
sample, i.e. for < x < L, the chemical potential is 
tuned to the Dirac point — 0). The function V{r) = 
^^^i Vn{r) represents N isolated scatterers. 

To calculate the full counting statistics of electron 
transport we use the matrix Green function approach 
[2^[2J,[23|. The Green function in the retarded-advanced 
(RA) space satisfies the equation 



where ^ = sin((/)/2) is the counting field and the chemical 
potential /i equals /Xqo in the leads and zero in the sample. 
It is important for the subsequent analysis that the bare 
Green function Qo, which solves Eq. ([!]) at V = 0, can be 
calculated analytically [23l |. 

Transport quantities (conductance, noise, and higher 
order cumulants) are readily determined from the cumu- 
lant generating function = Trln^""'^, where the trace 
Tr includes the spatial coordinates as well as RA, sublat- 
tice, and valley indices. Below we are mostly concerned 
with the conductance determined by the relation 



G = -{4eyh){d^T/d<P^) 



0=0' 



(2) 



With the help of the Dyson equation, Q^^ = Qq^ — V, we 
obtain J' = Fq+5T, where 5J^ = Tr \x\{l—V Qq) describes 
the correction to the cumulant generating function, J-q = 
Trln^o"^ = -Wcf)^ /2ttL, of the clean system ;23]. 

Our aim is to take advantage of the fact that impuri- 
ties do not overlap and to reduce the operator determi- 
nant in the definition of 5 J- to the usual matrix deter- 
minant. Such a reduction is justified if the impurity size 
is smaller than both the distance between impurities and 
the Fermi wave length (the latter is infinite at the Dirac 
point). The generating function can be rewritten in the 
TV-dimensional unfolded (impurity) space as 



5F =Tr\n{l-Vgo) = Tr ln(l - FSo), 



(3) 



where V = diag(Vi, V2, . . . , Vm) and all elements of the 
matrix Qq are identical and equal to the operator Q^. 

We proceed with introducing the T-matrix operators 
Tn = {I — Vng)^^Vn, that account for multiple rescat- 
tering on individual impurities 23:]. Here g stands for 
the Green function in an infinite graphene plane. At the 
Dirac point we find 

g(r,r') = -(i/27r) cr- (r - r')/|r - r'p. (4) 

Using the definition of the T-matrix, we rewrite SJ^ as 

SJ^^Tr\n[l-f{go-g)]+Tr\n{l-Vg), (5) 

where T = diag(ri, T2, . . . , Tjv) and g is proportional to 
the unit matrix in the unfolded space. 



The last term in Eq. ([5]) does not depend on the source 
field, (j), and can be safely omitted. Taking the limit of 
point-like impurity we reduce the operator product in the 
first term of Eq. ([5]) to the standard matrix product with 
the result 



^.F = Trln(l-f^,eg), 



(6) 



where T is the diagonal matrix of integrated impurity 
T-matrices (describing s-wave scattering only) and the 
elements of the matrix ^i-eg in the unfolded space are 
given by 



'f/o(r«,r,„). 



m ^ n, 



rcg J nm 



lim [t/o(r„, I") - .9(r„, r)] , m 



n. 



(7) 



where r„ = {xn,yn) specify the positions of impurities. 

Equation © is one of the central results of the Let- 
ter. The impurity-induced correction to the full-counting 
statistics is reduced to the finite-size matrix determinant, 
which is completely defined by the impurity T-matrices 
in the s-wave channel and the bare Green function. The 
T-matrices can be found in a standard way from the solu- 
tion of the corresponding single impurity scattering prob- 
lem [2^ . The exact bare Green function in the rectangu- 
lar geometry of Eq. ([T]) has been calculated in Ref. 23'] . 
It is given by the matrix product 

^o(r„, r,„) = ^t/,„AI],e^»(«"-«'")*/2L^^j;-i^ (g) 
where Tix.y.z are the Pauli matrices in RA space and 



(j){x — L) (p{x~L) 
= \ . Xx ■ ■ (hx 



A 



1 

(7, 



RA 



The matrix R is acting in the sublattice space. In the 
limit W ^ L it simplifies to 

f?(r r \ — / CSc(z,i + Z„) CSC(Z„ — Z„i)\ /„n 

\^csc(z„ - z„) csc(z„ + 2„0y ^ 

where z„ — 7r(a;„ -|- iyn)/2L and cscz — l/sinz. The 
matrix R coincides with the retarded Green function of 
the clean sample up to a factor AiL. 

The expressions ([8]), ^ define the off-diagonal ele- 
ments of the matrix Q^cg in the unfolded space. The 
diagonal elements are found from Eq. ([7]) as 



4i 



[/,„AE,i?,egAC/-i, 



(10) 



where i?i.og(r) — csc{ttx/L) + YiyOycji/iT. Diagonal part 
of the matrix i?rcg is proportional to the local density 
of states in a clean setup. Equation ([T0| completes the 
construction of the matrix 5rcg- 

The only input parameters for the general result (jS]) 
are T-matrices of individual impurities. They can be 
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obtained by solving the single impurity scattering prob- 
lem. For smooth potential impurity, the T-matrix mixes 
neither sublattices nor valleys and is given by the scat- 
tering length, T — £. The length £ diverges if the im- 
purity potential fulfills the resonant condition [S^l- In 
contrast, the T-matrix for an on-site potential impurity 
projects on a one-dimensional subspace, T = = 
£{1 =F TxCTx ± Tyay + TzOz)/'^, where upper (lower) sign 
correspond to the impurity in A (B) sublattice, respec- 
tively, and Ta are the Pauli matrices in the valley space 
[iij . The vacancy (i.e. infinitely strong on-site potential) 
corresponds to the limit £ ^ oo. 

In general, the T matrix of a resonant impurity is given 
by a divergent length scale, £, multiplied by a projection 
operator acting in the valley and sub-lattice space. This 
enables further simplification of the result ^ by omitting 
the unity under the logarithm in the projected basis. Up 
to an arbitrary constant term, the resulting generating 
function can be cast in the form 



6T = TT\TLkk\ 



(11) 



with a matrix K satisfying the identity K\4)) — k{—(j)). 
For resonant scalar impurities, the elements of 2N x 2N 
matrix K are given by 

fcr^i?(r„,r™), m/n, 

" i / /r^ ■ ^/ (1^^ 

\az csc{TTXn/ J^) — ^o■x4>/^T, m = n. 

For vacancies, the matrix K = A + A^ has a dimension 
N X N with 



A — 



exp 



l2L 



sin [ji;{CnXn + Cn 



0] 



(13) 



where = ±1 if the ji-th vacancy belongs to the sub- 
lattice A (B). The analytical expressions (fTT1) -(fT3 |) can 
now be used for the efficient numerical evaluation of the 
conductance, noise, and higher transmission cumulants 
of a disordered graphene sample. Below we focus on the 
conductance, Eq. ([2|). In terms of K it is given by 



4e2 r W 



27rTr 



{KK 



-1\2 



KR- 



4=0 



(14) 



where dots denote derivatives with respect to (j). In the 
case of resonant potential impurities, the matrix K is 
linear in (/>; hence the last term drops from Eq. (jl4p . 

Computational efficiency of Eq. (fT4| is limited by in- 
verting the matrix K &t (f) — {). This operation involves 
0{N^) multiplications. We run the standard matrix- 
inverse update algorithm by adding impurities one by 
one to compute the dependence of conductivity on N. 
This reduces the complexity to 0{N'^) per realization on 
average. The procedure is repeated many times to get 
sufficient statistics for different ratios W/L. Then we 
extrapolate the result to the limit — oo, thus elimi- 
nating non-universal boundary effects. 
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FIG. 1: (Color online) Mean conductivity calculated numer- 
ically from Eqs. (|12p . (|14p as a function of concentration of 
resonance scalar impurities (data points). The data is ob- 
tained by interpolation to the limit W ^ L. The solid line 
shows the result (|16p of the sigma model in class DHL 



For resonant potential impurities, the dependence of 
the average conductivity a = GL/W on the impurity 
concentration n = N /LW in the limit ly ^ L is plot- 
ted in Fig. [TJ To understand this behavior analytically, 
we perform the symmetry analysis of the problem. The 
matrix K defined by Eq. (|12p yields the Bogoliubov- 
de-Gennes type of symmetry, axK'^ {~(p)ax 



-K{q>), 

which corresponds to the symmetry class D 25|. The 



symmetry of the matrix K is equivalent to that of the 
transfer matrix of the system and can be used to infer the 
symmetry class of the corresponding Hamiltonian, which 
is given by Dili in the present case. The renormalization 
group analysis of the corresponding sigma model in the 
two-loop approximation yields the following equation for 
the dimensionless conductivity .26.] 



da 

d In L TT 



TTO" 



4e2 



(15) 



which holds for cr ^ 1, i.e., in the diffusive regime. Solv- 
ing Eq. (fT5|) . we get 



(Ae^/TTh) {\nnL^ - In In nL^) 



(16) 



Figure [T] shows a perfect agreement between the numer- 
ical and analytical results. 

An altogether different behavior of the conductivity is 
obtained from Eqs. (fT3|) . (fT4| for graphene with vacan- 
cies. The result is plotted in Fig.[2]in the limit W ^ L ioi 
different relative concentrations, ns/nA = 0, 1/2,3/4, 1, 
where ua and stand for the vacancy concentrations 
in the sublattice A and B, respectively. We see that the 
conductivity acquires a constant value for nL^ oo. 
To understand this behavior, we note that the matrix K 
from Eq. ([T3|) now possesses the only symmetry = K 
and thus belongs to the symmetry class Al. Therefore, 
graphene with randomly distributed vacancies falls into 
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FIG. 2: (Color online) Mean conductivity vs. vacancy concen- 
tration, n = JIA + nB, found numerically from Eqs. (|14|l . (I13p 
for ns/uA ~ 0,0.5,0.75, 1. Inset illustrates the conductivity 
scaling for ha ub on the logarithmic scale. 



the Hamiltonian symmetry class BDI. The corresponding 
sigma model is characterized by a vanishing /3-function, 
implying a constant (in general, nonuniversal) value of 
conductivity in the infrared limit. This is fully consis- 
tent with the numerical results of Fig. [2] Remarkably, 
for UA ^ ub, the conductivity is a non-monotonic func- 
tion of nL^ and the limiting value is very close to the 
4:6^ /nh. For equal concentrations, ha = ub, the numer- 
ically obtained limiting value is a* ~ 1.6 x Ae^/Trh. The 
system is expected to show various aspects of criticality 
characteristic for 2D problem of chiral classes [2^ . 

Additional comments: (i) If a small concentration of 
vacancies is added to the sample with resonant potential 
impurities, a crossover from Dili behavior, Eq. (|16p. to 
BDI (saturation) occurs at some high value of a. (u) We 
assumed that the inter-impurity distance is much larger 
than the graphene lattice constant a. For very large con- 
centration of potential impurities they will start to over- 
lap and will lose the resonant character. As to vacancies, 
when their concentration will increase towards ^ 1/a^, 
the conductivity will start to drop, and eventually the 
system will undergo a localization transition. 

In conclusion, we have developed a theoretical ap- 
proach to transport in disordered systems which de- 
scribes an entire crossover from ballistic to diffusive or 
critical regime. The theory can be applied to study local- 
ization physics and criticality in a variety of different sys- 
tems. We have used the theory to calculate the conduc- 
tivity (and, more generally, the full counting statistics) 
in undoped graphene with resonant impurities. The con- 
ductivity increases logarithmically in the case of smooth 
resonant potential scatterers (symmetry class Dili) and 
saturates at a constant value for vacancies (class BDI). 
In the latter case, the behavior of conductivity depends 
on the vacancy distribution among two sublattices. 
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[1] K. S.Novoselov et ai, Science 306, 666 (2004). 

[2] K. S.Novoselov et ai, Nature (London) 438, 197 (2005); 
Y.Zhang et al, Nature (London) 438, 201 (2005); ¥.- 
W.Tan et al, Eur. Phys. J. Spec. Top. 148, 15 (2007). 

[3] A. H. Castro Neto et al, Rev. Mod. Phys. 81, 109 (2009). 

[4] M. I. Katsnelson, Eur. Phys. J. B 51, 157 (2006). 

[5] J.Tworzydlo et al, Phys. Rev. Lett. 96, 246802 (2006); 
C. W. J. Beenakker, Rev. Mod. Phys. 80, 1337 (2008). 

[6] H. B. Heersche et al, Nature 446, 56 (2007). 

[7] F. Miao, et al. Science 317, 1530 (2007). 

[8] R.Danneau et al, Phys. Rev. Lett. 100, 196802 (2008). 

[9] K. I. Bolotin et al. Solid State Commun. 146, 351 (2008). 
[10] Xu Du, I. Skachko, and E. Y. Andrei, Int. J. Mod. Phys. 

B 22, 4579 (2008). 
[11] M. Titov, Europhys. Lett. 79, 17004 (2007). 
[12] J. H. Bardarson, et al, Phys. Rev. Lett. 99, 106801 
(2007); K.Nomura, M. Koshino, and S.Ryu, ihid. 99, 
146806 (2007); P. San- Jose, E. Prada, and D. S. Golubev, 
Phys. Rev. B 76, 195445 (2007); C. H. Lewenkopf, 
E. R. Mucciolo, and A. H. Castro Neto, ihid. 77, 
081410R (2008); J.Tworzydlo, C. W. Groth, and 

C. W. J. Beenakker, ibid. 78, 235438 (2008). 
A. Schuessler et al. Rev. B 79, 075405 (2009). 
Y.-W. Tan, Y. Zhang, H.L. Stormer, and P. Kim, Eur. 
Phys. J. Special Topics, 148, 15 (2007). 
P. M. Ostrovsky, I. V. Gornyi, and A.D.Mirlin, 
Rev. B 74, 235443 (2006). 

P. M. Ostrovsky, I. V. Gornyi, and A.D.Mirlin, 
Rev. Lett. 98, 256801 (2007). 
P. M. Ostrovsky, I. V. Gornyi, 
Phys. J. Spec. Top. 148, 63 (2007). 

T. O. Wehhng et al, Phys. Rev. B 75, 125425 (2007); 
ibid, B 80, 085428 (2009). 

T. Stauber, N.M.R. Peres, and F. Guinea, Phys. Rev. B 
76, 205423 (2007). 

D. C. Elias et al Science 323, 610 (2009). 
Z.H. Ni et al, arXiv:1003.0202 

J. Bardarson, M. Titov, and P. W. Brouwer, Phys. Rev. 
Lett. 102, 226803 (2009). 

M. Titov et al, Phys. Rev. Lett. 104, 076802 (2010). 
S.Ryu et al, Phys. Rev. B 75, 205344 (2007). 
for review see F. Evers and A.D. Mirlin, Rev. Mod. Phys. 
80, 1355 (2008). 

In fact, the sigma model contains an additional Wess- 
Zumino term (in analogy with Z2 topological term in 
class All for non-resonant potential scatterers [3). How- 
ever, it does not affect the two- loop /^-function at large 
a (I15|) and thus the asymptotic behavior of conductivity 
(|16|l . The /3-function has a zero ai a — I/tt corresponding 
to the Wess-Zumino-Witten fixed point. A remarkable 
property of class DIII is that this fixed point is unstable: 
a larger than I/tt flows to infinity. 
[27] Yu. V. Nazarov, Phys. Rev. Lett. 73, 134 (1994). 
[28] M. Hentschel and F. Guinea, Phys. Rev. B 76, 115407 
(2007); D. S. Novikov, Phys. Rev. B 76, 245435 (2007); 
D.M.Basko, Phys. Rev. B 78, 115432 (2008). 



[13 
[14 

[15: 

[16 

[17 

[18 

[19 

[20 
[21 
[22; 

[23' 
[24' 
[25 

[26 



Phys. 
Phys. 

and A.D.Mirlin, Eur. 



